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NOTATION 


ROMAN 


1 . Italicized Lower Case 


9 

k 

1 

m 

n 

9 

r 

s 

t 

V 

X 


scalar mass flux 
thermal conductivity 
length 
index limit 
index limit 

contact heat flux scalar 
radius 

condensation 

time 

velocity scalar 
displacement scalar 


2 . Italicized Upper Case 


A 

E 

K 

M 

N 

P 

Q 

R 

T 

U 

V 


area 

external and mutual energy 

constant 

mass 

non-dimensional parameter 

pressure 

heat 

gas constant 
temperature 
internal energy 
volume 


3 . Bold Italicized Lover Case 
g mass flux density 


ii 



n 

unit outvard 

normal 

q 

contact heat 

flux 

V 

velocity 


X 

displacement 



4 . Bold Upper Case 

T extra stress tensor 

5 . Lover Case 

g gravitational acceleration 

p particle 

GREEK 

1 . Upper Case 

E entropy 

T mole fraction 

2 . Lover Case 

a constant 

k thermal diffusivity 

H dynamic viscosity 

v kinematic viscosity 

p density 

generalized scalar, vector, or tensor quantity 
cj angular velocity 


OPERATORS 

d total derivative 


iii 


D 


f() 


a 

A 

V 

J 

£ 


sub s t ant i ve de r i va t ive 

function of 

partial derivative 

incremental change 

divergence 

integral 

summation 


* it]i> time average of 
{V ]^ volume average of 

time average of volume average of 
j^l absolute value or magnitude of ^ 

* scalar product of vectors, vector product of vector and tensor 

: scalar product of tensors 


SUBSCRIPTS 


a 

b 

c 

ch 

comp 

eff 

h 

1 

is 

J 

k 

K 

Ha 

(m) 

n 

Pr 

Re 

(s) 

sys 

V 

uw 


acoustic 

boundary 

cold cavity 

characteristic 

computational 

effective 

hot cavity 

index 

isolated 

index 

index 

Kurzweg 

Mach 

material body 

momentum discrete volume 

Prandtl 

Reynolds 

system of particles 
system 

at constant volume 
upwind 


IV 



Va 

0 


Valensi 

fiduciary 

at constant entropy 


SUPERSCRIPTS 


s previous time step 

* distinguishing indicator 
per unit mass 

* time rate of change 

j index 



INTRODUCTION 


The activities described in this report cover two aspects pertinent to 
Stirling machine simulation, namely, information propagation and entropy 
transport. Information propagation issues have been shown in previous work to 
manifest themselves in hardware such as transmission lines (Go90) . While such 
effects have not been definitely demonstrated for Stirling machines, there are 
enough suggestions (Go87 , Go90) that such effects may be of concern in high 
frequency, high working pressure devices such as the Space Power Demonstrator 
Engine (SPDE) to warrant further study. Hence the work covered in this report 
is aimed essentially at attempting to determine whether the discretised 
primitive (or first order temporal differential) conservation balances are 
adequate to simulate information propagation effects in circumstances when 
they are known to be of importance. 

Simulation of entropy in Stirling machines offers the opportunity of 
giving the designer new insights into the sources of irreversibilities as well 
as a criterion for judging the effectiveness of design optimizations. 

Previous oscillating flow work in this area has been effective in simulating 
the generation of entropy (Ko90) while less attention has been focussed on its 
transport through the computational domain. The work reported herein thus is 
aimed at correcting this imbalance by focussing on the transport of entropy 
which includes a means of constructing the computation domain as an isolated 
thermodynamic system. 

Rather than simulating Stirling hardware directly, the approach adopted 
has been to base the simulation investigations on simple physical systems for 
which experimentally validated closed- form analytic solutions exist. 
Specifically, the information propagation study is based upon Iberall's 
analytic solution of the transmissi on li ne problem (Ib50) while the entropy 
transport investigation is carried out in terms of Kurzweg's oscillating flow 
apparatus (KZ84) . This approach enables the behavior of the simulation in 
terms of information propagation and_ entropy transport to be studied in 
isolation. This avoids the complicating gas-dynamic interactions inherent In 
Stirling hardware from masking the effects under study and consequently 
leading to erroneous conclusions. 

The report is divided into two parts. Part A describes the information 
propagation investigations while part B is devoted to the entropy transport 
work. 
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PART A 


INFORMATION PROPAGATION 


A. 1 INTRODUCTION 


The investigation of information propagation effects is carried out in 
terms of the following three tasks, namely: 

a. The development of a new, non- pressure linked numerical algorithm for 
applying the mass, momentum and energy transport equations in one 
dimension and the application of the algorithm to Iberall's transmission 
line problem. 

b . An investigation of the use of a hybrid explicit/ implicit integration 
algorithm in Stirling machine simulation using the transmission line as 
a validation vehicle. 

c. The development of an n-step (Go90) implicit temporal integration 
algorithm and its application to the transmission line to test whether 
the numerical integration algorithm itself influences the simulated 
Information propagation behavior. 

The analytic and physical details of Iberall's transmission line are 
described in detail in a previous NASA contractor report (no, 185285) (Go90) 
and hence will not be repeated here. 

A. 2 ANALYTIC BACKGROUND 


The essence of the pressure -linked algorithm is obtained by substituting 
the momentum equation into the continuity equation In which an equation of 
state is used to express density in terms of pressure. Following the 
derivation presented in Go90, the resulting equation may be expressed as: 


[tv]*? 
R [tvl TAt 


At 


72 (S) 


f- f- [tV] PndA • n 

J A„ ( „ l J A (j) J 
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•ro'Jl (l»‘ ■ -■)" 


(A. 2. 1.1) 


where : 
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(A. 2. 1.2) 


([tv n ]g v n(S)) S + At f Qtv„]g> [tv] r ' [tv]/’ ) 

V THS) 


This equation yields an advanced time or implicit pressure field. The 
information propagation characteristics of equation (A. 2.1) may be 
demonstrated by applying it to a sequence of adjacent discrete volumes in a 
one -dimensional field of constant cross sectional area. Hence simplifying the 
equation and dropping the averaging notation for the sake of clarity yields: 


ii 


A V" &t 2 A 2 RT (j> . . 
+ L, ^2 


1J 


) = K 


source 


(A. 2. 2) 


where is the diagonal or "central" pressure term, P±j are the outlying 

pressures and ^ sour ce the source term. Now V/A represents the length of 
the discrete volume in the computational direction and hence V/AAt represents 
a computational speed v comp since At is the integration time increment. Thus 
equation (A. 2. 2) may be expressed as: 


E 
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r comp 


E 
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^2 


f comp 


±3 


v source 


(A. 2.3) 


where v a is the isothermal speed of sound for a fluid with constant specific 
heats. Thus, in simple terms, the information propagation is determined by 
the ratio v a /v comp (acoustic speed to computational speed) . This conveniently 
explains the behavior of the pressure linked analysis as it applies to the 
transmission line as well as, for that matter, to the SPDE (Go90, chapters 3 
and 4). In particular, for the transmission line, it was observed that the 
transducer cavity to excitation pressure ratio asymptoted towards the value 
calculated via Iberall 's analysis as the time step was decreased, which is a 
classic numerical result. Simultaneously, the phase angle between the 
excitation and transducer cavity started out at a value much higher than the 
Iberall prediction and decreased towards an asymptote with a value 
approximately half that predicted by Iberall as the integration time step was 
decreased. This phase angle behavior is explained by equation (A. 2. 3) since 
simply increasing v comp by decreasing At decreases the coupling between p jj 
and P_£_£ until the Courant limit defined by v a - v C omp is reac ^ed. Under 
these conditions, the phase angle predicted by equation (A. 2.1) should 
correspond with the asymptotic value produced by a simulation in which the 
pressure field is extracted explicitly from the continuity equation, that is, 
without the use of any information propagation (or implicit pressure field) 
equation at all. Under these conditions, the requirement for using a unitary 
pressure domain (UPD) or equilibrium hypothesis for modelling information 
propagation (see Go90 section 2.6) would fall away and the criterion for 
determining the time step size is the achievement of numerical stability. 

Noting that the information propagation modelling capability per se is 
quantified by the phase angle predictions, the above observations lead to the 
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proposition that the discretised primitive conservation balances intrinsically 
are deficient as a basis for simulating the information propagation effects 
occurring in a transmission line (specifically , low Mach number or stationary 
flow resonance induced shocks (JI73)). Furthermore, by extension, this also 
mav be true for the information propagation effects occurring in Stirling 
machines with characteristic numbers less than 24. These effects include the 
superposition of multiple reflected waves and microscopic scale regenerator 
choking. 

The principal argument in support of the proposition arises from 
Iberall's analysis itself since Iberall does not solve the primitive 
conservation balances but in fact solves the Kirchoff equations of sound 
(Ib50, page 100, section 4). In particular, the relevant information 
propagation equation is obtained by substituting a simplified differential 
mass balance into a differential momentum balance which yields (quoting 
Rayleigh (St45) article 348, equation (2)): 

dv + _1_ dP = Jf_v 2 v + A* d2 f* (A. 2 . 4) 

"cTt p 0 "3x p 0 3 p 0 dxdt 

where s Is the "condensation" defined by: 

s = P SJl (A. 2. 5) 

P 0 

and p 0 is a constant, fiduciary density. 

Equation (A. 2. 4) is the inverse symmetrical form of the well-known wave 
equation obtained by substituting the differential momentum balance into the 
differential mass balance. As noted by Organ (OJ89) , a simplified isothermal 
form of the wave equation is given by: 

_Lif£=V z P (A. 2. 6) 

v 

which has the characteristics solution: 

P = £(t ± JL) (A. 2. 7) 

v a 

Hence in respect of accurately modelling information propagation 
effects, Iberall in effect agrees with Organ that the equations of sound 
should be used to extract the pressure fie ld and not the primitive 
conservation balances. In the past, Organ has attempted. tq solve these 
equations using the method of characteristics (Or82) , but, recognizing the _ _ 
Intractability of this approach for the geometrical complexities of Stirling 
machines, he has resorted to a linear solution of equation (A. 2. 6) using the 
principle of superposition (OJ89). Implicit In this argument is the notion 
that the second order differential equation (A. 2. 6) is not equivalent (at 
least in numerically discretised terms) to the first order pressure-linked 
equation (A. 2.1) even though their method of derivation is ostensibly the same 
(that is, substitution of a momentum equation into a continuity equation). 
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Rigorously, this can be shown to be the case since the integral transform of 
the complete version of equation (A, 2.6) obtained via the generalized 
transport equation (Go90, equation (A. 2. 3)) is very different from equation 
(A. 2.1). Therefore, this argument offers a reasonable explanation for the 
discrepancies observed to date between the predictions made by the simulation 
using equation (A. 2.1) and Iberall's analysis - the simulation is based upon 
an effectively different description of information propagation. 

However, if it is indeed true that the discretised primitive 
conservation balances used in the simulation inherently are incapable of 
describing the information propagation effects under discussion in this 
report, then all possible solutions of the primitive equations should show the 
same behavior, not just the pressure -linked forms evaluated to date. In 
particular, the symmetry between the second order differential Kirchoff 
momentum equation and the second order differential wave equation discussed 
above must also hold for the first or der integral equations. Thus a solution 
obtained from the symmetrical form of the pressure-linked equation (A. 2.1) not 
involving an implicit solution of the pressure field should show the same 
behavior as the pressure-linked form. 

Therefore, if both symmetry and explicit solution bounding of an 
implicit solution can be demonstrated for the transmission line, then the 
proposition that the primitive discretised conservation balances are 
inadequate for modelling transmission line information propagation effects 
would gain significant substance. 

A third element involved In this process is the notion that the 
qualitative information propagation modelling behavior of the primitive 
equations should be numerical integration process independent, particularly 
for the implicit forms. (It can be taken as read that this independence has 
been shown for explicit numerical integration, for example, see Roach (Ro82)). 
A means of testing this assertion is to implement the n-step implicit temporal 
integration scheme outlined in Go90 and apply it in both globally implicit and 
conventional (a matrix inversion at each time step) environments. In both 
cases, therefore, if the proposition is true, then the time step dependent 
phase angle behavior observed to date should be invariant. 

The symmetry test is conducted in task a, the explicit solution bounding 
test forms task b while the implicit numerical invariancy test is carried out 
as task c. 


A. 3 A NON -PRES SURE -LINKED INVERSE SYMMETRIC INTEGRAL FORMULATION 

If the pressure-linked formulation of equation (A. 2.1) is obtained by 
substituting the momentum conservation balance into the continuity equation, 
then the inverse symmetric form is obtained by substituting the continuity 
equation into the momentum conservation balance, so eliminating the pressures. 
Hence substituting the equation of state (Go90, equation (2.33)) into the 
continuity equation (Go90, equation (2.41)), temporally discretising the left 
hand side and rearranging produces: 
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Substituting equation (A. 3.1) into the momentum equation (Go90, equation 
(2.42)) yields: 
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(A. 3.2) 


where g * represents the laminar and turbulent shear stress terms. Equation 

(A. 3. 2) yields the advanced time or implicit mass flux field g. As before, 
the information propagation characteristics of the equation may be 
demonstrated by considering a sequence of adjacent discrete volumes in a one- 
dimensional field of constant cross -sectional area. Thus, temporally 
discretising the left hand side of equation (A. 3. 2), simplifying and dropping 
the averaging notation yields: 

+ a £{( gv) u - (gv)ij) * £ £*™± (g u - gij ) - K soujrce . 


Rearranging, and noting that v ( 


comp 


V/AAt and v 2 - RT yields: 


igv) 
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r ~ ^ 

- 

r 

v comp + ^ 

i . 1 v ?„ 

- y 

i * 1 v * 


**Ma v comp 

J j 

^ Ma ^ comp 


(gv) 


j-j 


(A. 3. 3) 


v source 


where is the Mach number. Therefore for the transmission line, the 

Vo 


dominant term is 


v comp N Ma 


since N Ma « 1. This is also true in the Courant 


limit when v a = v C omp’ Hence the " information propagation '* characteristics 
(in terms of mass flux) of equation (A. 3.3) have a discretised description 
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different from those of the pressure -linked equation (A. 2.1) and thus use of 
equation (A. 3. 2) provides a valid test for evaluating the symmetry of the 
primitive conservation balances. 

Thus the set of steps comprising the solution algorithm becomes: 


a. Guess the mass flux field [tVn3 g* . 

jb. Explicitly compute the discrete volume masses from Go90 equation (2.41) 
and then infer [tV] p from known values of V {3) . 

c . Implicitly compute the temperature field using Go90 equation (2.43). 

d. Explicitly compute the pressures from the equation of state (Go90 
equation (2.33)). 

e. Implicitly compute the mass flux field from equation (A, 3. 2). 

f. Compare [tVn ]9 (computed) with [t v n ]9* (guessed) and return to step b 

with f ( [tv n ]9> [tv n i 9* if the mass flux field is insufficiently 

converged. 

A. 4 EXPLICIT SOLUTION BOUNDING FORM 

The hybrid explicit/implicit algorithm proposed by MacCormack (Mc82) 
originally was Investigated with a view to improving the efficiency of the 
existing purely implicit algorithm since it c^fensibly offers a means of 
eliminating the iteration procedure on which the algorithm is based. 
MacCormack' s algorithm may be viewed as a combination of two Independent 
entities. The first is the hybrid explicit/implicit integration algorithm 
itself, while the second is an elaborate factorization procedure for 
extracting a pseudo -characteristics solution from the Navier- Stokes equations 
without explicitly solving an M Information propagation 1 ' equation of the 
pressure- linked variety. However, the factorization procedure does not appear 
amenable to ready generalization (at least in discrete volume, integral terms) 
since it is derived for a conformal mesh (Cartesian in Mc82) using the partial 
differential form of the conservation equations. 

In view of the discussion given in section A. 2 above, only the hybrid 
explicit/implicit numerical algorithm portion of MacCormack' s algorithm has 
been investigated since the pseudo-characteristics portion preempts an 
Isolated evaluation of the primitive conservation balance inadequacy 
proposition. In other words, if MacCormack' s full algorithm were to replicate 
Iberall's analytic results, it would not be possible to determine whether this 
was due to the implicit/explicit algorithm or the pseudo-characteristics 


7 


ORIGINAL PAGE IS 
OF POOR QUALITY 


formulation. Thus in view of time and funding constraints, investigation of 
this latter aspect of MacCormack's algorithm as applied to a transmission line 
was not attempted. 

The basic elements of MacCormack's hybrid explicit/implicit numerical 
integration algorithm (as applied to the integral conservation balances) may 
be described by the following sequence in terms of a general variable rf>: 

a. Explicitly calculate the temporal gradient at the current time step: 

**(t) - f(t f *(t)) (A.A.l) 

jb. Implicitly extract an estimated advanced time value from: 

i|jf!_(t+At) = (l-a)i^(t) + a f (t + At , ( t + At) ) (A. 4. 2) 

dt dt 

c. Explicitly calculate the estimated advanced time temporal gradient: 

( t + At ) = f (t + At , ( t + At) ) (A. 4. 3) 

dt 

d. Implicitly extract a corrected advanced time value ^ from; 

^(t+At) = (1 - a) (t + At) + a f (t + At , if>(t +At) ) (A. 4. 4) 

e. Assemble an averaged final advanced time value: 

V>(t + At) = + -^At ^ (t +At) + ||(t + At)j (A. 4. 5) 

♦ 

In steps b and c, a is referred to as an "implicit blending factor". 

On closer examination, the algorithm can be seen to be an interesting 
assembly of several existing methods. When a - . 5 (the numerical stability 
limit), steps a and b comprise a single iteration of the Crank-Nicholson 
algorithm while steps c and d comprise a single Iteration of the standard 
implicit algorithm (in this case, a corresponds to the convergence factor). 

Step e represents an aggregation procedure similar to that found in the Runge- 
Kutta algorithm and its derivatives (such as those developed by Fehlberg and 
Verner) , 

The algorithm was applied to the standard set of pressure -linked 
Integral equations and applied to the transmission line using the same 
parameters listed in Go90, section 4.9. In the first instance, the 
convergence of the algorithm at relatively large time steps (such as those 
generated by the UPD hypothesis) was quite poor (that is, large truncation 
errors). Hence several Iterations of steps c, c? and e were necessary in order 
to produce a stable solution. These Iterations submerged the benefit of steps 
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a and b t In effect making them redundant. Furthermore, there was no change In 
the observed phase angle or amplitude ratio behavior in comparison with that 
shown by the standard Crank-Nicholson algorithm (Go90, sections 4.8 and 4.9). 

Hence, in isolation, the particular explicit/ implicit formulation of 
MacCormack does not contribute towards an understanding of the information 
propagation issues under discussion, in particular, it does not change the 
influence of the pressure -linked equation (A. 2.1). Furthermore, it does not 
appear to yield any computational advantage in comparison with the standard 
iterative implicit method. However, it must be emphasized that these 
observations may not be true for MacCormack* s method as a whole, that is, 
including the pseudo -characteristics formulation of the conservation 
equations . 

Therefore, in order to ascertain the limiting behavior of the primitive 
conservation balances without the use of a pressure- linked formulation (as 
discussed in section 2), the standard implicit algorithm (Go90, section 2.5) 
is modified to remove the pressure- linked equation (A. 2.1) so that the 
pressure field is extracted via an equation of state from the explicit mass 
conservation balance. The modified algorithm may be listed as follows: 


a. Guess the mass flux field [ t v n )9* * 

j b. Explicitly compute the discrete volume masses from Go90 equation (2.41) 

and then Infer [tV] p from known values of V (S) . 

c. Explicitly compute the pressures from the equation of state (Go90 
equation (2.33)) using current estimates of the advanced time 
temperature field. 

d. Implicitly compute the temperature field using Go90 equation (2.43). 

e. Implicitly compute the mass flux field from (Go90) equation (2.42) using 
the explicitly derived pressures of step c as part of the source term. 

f . Compare [tv n }9 (computed) with [tVnj g* (guessed) and return to step b 

with f( [tVn] g, {t v n ]9*) lt v n j 9* if the mass flux fie l d is Insufficiently 

converged . 

This algorithm allows the explicit information propagation boundary of 
the implicitly solved primitive conservation balances to be determined without 
introducing any complicating numerical effects. 
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A. 5 n - STEP IMPLICIT TEMPORAL INTEGRATION ALGORITHM 


The n - step Implicit algorithm is based upon a simple multi- increment 
expansion of the left hand side of the integral conservation balances. 
Consider a generalized form for the transport property rjr. 


d (tV] tfV 


[tvjtf ) 


(A. 5.1) 


Discretising the left hand side linearly in terms of n time steps yields: 


(ltvi^ v ) J " ([tv]^^) J 1 = f 
~KtJh 


+ 2^1 


n 


» [tV] 


iP 


(A. 5. 2) 


where j is the time step index. Clearly, a linear expansion is not a 
necessary restriction, other expansions using, for example, sinusoidal or 
polynomial forms are equally valid. However, use of more complicated 
expansions does not affect the argument in principle, although such higher- 
order expansions probably yield superior numerical results. 


Hence, in one -dimensional terms, if there are m discrete volumes, the 

solution is obtained from a system of mxn equations, each with the implicit 
form (after dropping the averaging notation for clarity) : 


• V) ± 


(A. 5. 3) 


where i is the spatial discretisation index. For the first time increment 
(j — 1), eqiaat: ion (A. 5, 3) becomes: 


■<**)£, + (**){ - « (*i) 


source 


A t/n 


(A. 5.4) 


where the superscript s denotes the initial conditions. 

In general, there Is no restriction on At , m or n so that At may be 
equated with the cycle period. Under these conditions, the algorithm becomes 
analogous to Gedeon's "globally implicit" technique with the important 

exception of equation (A. 5. 4), that is, the assumption that = not 

made. This version of the globally implicit technique yields a transient 
evolution of the flow field without imposing external cyclic closure boundary 
conditions. Furthermore L J:he n - step formula tion is independent of the 
specific form of equation (A. 5. 1) . Hence It can be applied to the pressure- 
linked equation (A. 2.1), the symmetric momentum linked equation (A. 3.2) as 
well as to the unmodified momentum and energy conservation balances (Go90, 
equations (2.42) and (2.43) respectively). Thus the iterative integration 
algorithms described in sections A. 3 and A. 4 and in Go90 section 2.5 may be 
used without alteration using the n - step formulation. This allows the 
integration algorithm independence of the primitive equations to be assessed 
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without adulteration by extraneous effects as required by the discussion in 
section A , 2 . 

A. 6 RESULTS 

The simulations are tested for two transmission line geometries and 
their respective boundary and excitation conditions as detailed in table A.l 
giving a total of five test cases. The geometries and boundary conditions 
are identical to two of the test cases reported by Watts (Wa65) . Particular 
care was taken to use exactly the same excitation pressure amplitudes as used 
experimentally. Where possible, the excitation frequencies were chosen to 
correspond with experimental tests in order to provide an experimental % 

baseline for the comparison of the simulation and experimental data. In this 
regard, it must be pointed out that neither Watts nor Goldschmeid (G 068 ) 
report any experimental data for validating the phase angle predictions of 
equation (116) in Ib50. Hence it is assumed here that because Iberall's 
amplitude ratio predictions (also obtained from equation (116)) are confirmed 
by experiment, his phase angle predictions are likewise valid. However, this, 
may not be true , a qualification which must be borne in mind when reviewing 
the results presented below. 

The six inch line excited at 100 Hz has a characteristic number 1 (tf c73 ) 
of 22.86 (see table A. 2) and is equivalent in information propagation terms to 
the SPRE which has N ch - 24. Hence to the extent that the information 
propagation effects occurring in the transmission line match those occurring 
in the SPRE (as alluded to In section A, 2, this similarity may be only 
nominal), the simulation comparisons provide a test of the information 
propagation modelling validity of the simulations as applied to the SPRE. 

The results for the five test cases are given in tables A. 2 through A , 6 
located at the end of this section and are presented in terms of amplitude 
ratio and phase angle. The amplitude ratio Is defined as the ratio between 
the transducer cavity and excitation pressure amplitudes while the phase angle 
is the phase lag between the excitation and transducer cavity waveforms. With 
reference to tables A. 2 through A. 6 , the first line reports the analytic 
predictions of Iberall's analysis and the second line gives the amplitude 
ratio measured by Watts (Wa65) (where available). Thereafter, the tabulation 
gives the results of eight different simulations each corresponding to a 
different integration algorithm / equation formulation mix. 

The first and third simulations correspond to the UPD and equilibrium 
information propagation hypotheses (defined in Go90, section 2.6) using the 
pressure - linked equation (A. 2.1) (obtained by solving the continuity equation 
after substitution of the momentum equation, hence the "UPD: continuity" and 
"Equilibrium: continuity" nomenclature). 


1 The characteristic number N ch is defined as the number of pressure wave 
traverses occurring per cycle over the entire extent of the working space. 
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Table A.l Transmission line geometries 


Parameter 

6 inch line 

24 inch line 

Length (mm) 

152.4 

609.6 

Diameter (mm) 

4.66 

4.66 

Cavity volume (cm 3 ) 

.414 

.414 

Gas type 

air 

air 

Mean pressure (bar) 

.77 

.77 

Excitation pressure 
amplitude (bar) 

.00077 

.00077 

Temperature (°C) 

28.9814 

25.95 

Excitation frequencies 
(Hz) 

100; 230; 400 

55; 110 


The second and fourth simulations portray the UPD and equilibrium information 
propagation hypotheses using the symmetric momentum equation (A. 3. 2) (obtained 
by solving the momentum equation after substitution of the continuity equation 
yielding the "UPD: momentum" and "Equilibrium: momentum" nomenclature). Both 
equilibrium simulations are performed at 50 and 100 increments per cycle in 
order to confirm the phase angle and amplitude ratio behavior. The first four 
simulations are thus designed to accomplish the symmetry test. 

The fifth simulation uses the algorithm described in section A. 4 
(explicit solution bounding) in which the momentum and energy balances are 
solved implicitly and the pressure field is determined explicitly without a 
pressure -linked equation. The number of increments per cycle is chosen purely 
on the basis of numerical stability, being small enough to achieve a stable 
solution but not so small as to undercut the Courant criterion and introduce 
dispersion errors (Ro82). 

The last three simulations test various forms of the n - step implicit 
algorithm. The sixth simulation is based on the UPD hypothesis for 
determining the number of time steps per cycle, while each time step is 
divided into n substeps. Four values of n are tested in order to establish 
the behavioral trends of the amplitude ratio and phase angle. The seventh 
simulation embodies the globally implicit version of the n - step implicit 
algorithm while the last simulation is similar to the "UPD: n - step implicit" 
simulation except that the time steps are variable and are determined 
characteristically at each stage of the simulation. 

The results for the test case approximating the SPRE in transmission 
line information propagation terms are given in table A. 2 (6 inch length with 
100 Hz excitation) . This case yields an analytic characteristic number of 
22.862 which is confirmed within a tolerance of +.002 by all the simulations. 


12 


























Iberall's analysis predicts an amplitude ratio close to unity and .03 higher 
than that measured (the measured value is read manually off a logarithmic 
graph in Wa65 and therefore is only approximate) while the predicted phase 
angle of .5° is small. The predictions made by the momentum and continuity 
linked equations are identical and show the same trends as the time step size 
is decreased (number of increments per cycle is increased). As noted before, 
these trends are typified by the amplitude ratio increasing towards the 
Iberall value as the time step size is decreased while, concomitantly, the 
phase angle decreases. Hence, in this case, the symmetry of the equations is 
demonstrated. 

The simulation with no pressure linking using 300 increments per cycle 
produces an amplitude ratio equal to that produced by the 100 increment per 
cycle equilibrium algorithm (at least to three decimal places) and a phase 
angle lower than that produced at 100 increments per cycle. Hence, in this 
case as well, an explicit determination of pressure indeed does produce a 
bounding value of amplitude ratio and phase angle for the implicit continuity 
and momentum linked equation algorithms. 

Without going into needless repetitive detail, these observations apply 
for the remaining test cases as well (tables A. 3 through A. 6). A few 
interesting points to note are: 

In table A. 4, at 400 Hz excitation, Iberall's analysis underpredicts the 
experimental amplitude ratio by 31% while the asymptotic simulated 
pressure amplitude ratio for the first five simulations of 3.301 is 15% 
less than Iberall's value and 41% less than that measured. 

Also in table A. 4, the equilibrium simulations as well as the case with 
no pressure linking show identical results at a temporal discretisation 
of 100 increments per cycle. This further demonstrates that the 
explicit pressure determination defi nes the bounding case of the 
implicit determination using either continuity or momentum equation 
linking. 

In table A, 6, a comparison of the predictions of the equilibrium 
simulations with the predictions produced by the simulation with no 
pressure linking shows that when the temporal discretisation of the 
explicit pressure simulation is coarser than that of the implicit 
simulations, then the implicit simulations predict phase angles smaller 
and amplitude ratios larger than the explicit simulation. This conforms 
to the standard time step size dependent trend and shows that at small 
enough time steps, extracting the pressure field implicitly via the 
continuity and momentum linked equations is computationally redundant. 

Therefore, the first five simulations demonstrate that the continuity 
and momentum equation linked equations are indeed symmetric and that the 
implicitly extracted pressure fields are the same as those extracted 
explicitly at small enough time steps. In other words, explicit pressure 
bounding for the primitive conservation balances Is confirmed. 
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Consequently, the remaining n step implicit sim ulations were carried 
out using the momentum linked formulation as it is numerically more efficient 
than its symmetrical continuity or pressure - linked counterpart because one 
less matrix inversion is required for each iteration. 

Reverting to table A. 2, it may be observed that the UPD and 
characteristically determined time step versions (numbers 6 and 8) produce 
identical results in terms of amplitude ratio versus increments per time step. 
The trend in both cases conforms to the standard pattern, namely, amplitude 
ratio increasing with decreasing sub- time step size. Furthermore, the 
amplitude ratio at 13 increments per time steps conforms to the limiting value 
produced by the single step implicit simulations (numbers 3 and 4) . In terras 
of phase angle behavior, the UPD and characteristic versions of the n - step 
implicit algorithm show the same trend (phase angle decreasing with increasing 
number of increments) except that there is some disparity between the 
numerical values. This is caused by the relative lack of cyclic convergence 
of the characteristic version (cyclic energy balance errors > 1 % versus 
equivalent errors for the UPD version < .01%). The characteristic version 
takes an inordinately large number of cycles to reach cyclic equilibrium and 
hence the simulation was te rminat ed before full convergence was achieved as a 
matter of expediency. Hence the characteris tic version, while being a useful 
test of numerical algorithm Independence, is not a practical methodology. 

Finally, the behavior of the n - step globally Implicit algorithm shown 
in table A, 2 is rather interesting. A discretisation of 20 increments per 
cycle proved to be the largest practical value that would enable the 
coefficient matrices to be Inverted (using a computer with a 64 bit precision 
word length) before numerical round-off errors caused the simulation to become 
unstable. This situation could be alleviated somewhat by various 
normalization schemes as well as by matrix partitioning algorithms. However, 
In the context of this investigation, such endeavors would not add materially 
to the outcome and hence were not attempted. In this context, the n - step 
globally implicit scheme conforms to the established pattern of amplitude 
ratio and phase angle behavior observed for all the other simulations. 
Examining the phase angle behavior, one is tempted to infer that the value 
asymptotes to the Iberall value of .525°. However, in temporal discretisation 
terms, at 20 increments per cycle, the algorithm is In the range of the UPD 
algorithm, which with = 22.864 yields a temporal discretisation of 23 
increments per cycle. Thus at 20 Increments per cycle, the n - step globally 
implicit algorithm produces a phase angle greater than the 23 increment per 
cycle UPD algorithm which is In conformity with the established trend. 

With minor variations, these observations also apply to the results 
produced by the n - step implicit simulations for all the remaining test cases 
in tables A. 3 through A. 6. The only exceptions worth noting occur in tables 
A. 4 and A. 6 in which the phase angle produced by the globally Implicit 
simulation at 10 increments per cycle is greater than that produced at 5 
increments per cycle. Once again, these exceptions are caused by cyclic 
equilibrium convergence effects where the particular simulations are not as 
converged as their counterparts. However, since the phase angles in all cases 
at 20 increments per cycle are less than those at coarser discretisations, the 
basic behavioral trend is not violated. 
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Therefore, the n - step globally implicit simulations confirm that the 
simulation predictions are independent of the particular form of the 
integration algorithm used to solve the conservation balances. Furthermore, 
use of an n - step implicit algorithm does not appear to offer any numerical 
accuracy advantages that cannot be achieved at a lower computational cost 
using a single step implicit Crank-Nicholson algorithm with a finer temporal 
discretisation. In this vein, there does not appear to be any advantage in 
using a globally implicit algorithm, at least in terms of a strict transient 
simulation. 
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Table A. 2 6 Inch transmission line with 100 Hz excitation 


Reference 

Number 

Analysis 

Temporal 

Discretisation 

Parameter 

Parameter 

Value 

Char . 
Number 

Amplitude 

Ratio 

Phase 

Angle 

(deg) 


Iberall 

(analytic) 

-- 

— 

22.862 

1.060 

0.525 


Watts 
(exper . ) 


-- 


1.03 

-- 

i 

UPD: 

continuity 

-- 

-- 

22.864 

1.049 

0.97 

2 

UPD: 

momentum 

-- 

-- 

22.864 

1.049 

0.97 

3 

Equilibrium: 

continuity 

Increments 

per 

cycle 

50 

22.864 

1.051 

0.59 

100 

22.864 

1.052 

0.41 

4 

Equilibrium: 

momentum 

Increments 

per 

cycle 

50 

22.864 

1.051 

0.59 

100 

22.864 

1.052 

0.41 

5 

Implicit : 
no pressure 
linking 

Increments 

per 

cycle 

300 

22.864 

1.052 

0.32 

6 

UPD: 

n - step 
implicit 

i 

Increments 

per 

time step 

2 

22.864 

1.050 

0,80 

3 

22.864 

1.051 

0.56 

5 

22.864 

1.052 

0.38 

13 

22.864 

1.052 

0.31 

7 

n - step 
globally 
implicit 

Increments 

per 

cycle 

5 

22.864 

1.011 j 

2.69 

10 

22.864 

1.038 

1.98 

20 

22.864 

1.048 

1.11 

8 

n - step 
implicit : 
character- 
istically 
determined 
time step 


2 

22.860 

1.050 

0.54 

Increments 

3 

22.860 

1.051 

0.54 

per 

time step 

5 


1.052 

0.43 


13 

22.860 

1.052 

0.18 
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Table A. 3 6 inch transmission line with 230 Hz excitation 


Reference 

Number 

Analysis 

Temporal 

Discretisation 

Parameter 

Parameter 

Value 

Char. 

Number 

Amplitude 

Ratio 

Phase 

Angle 

(deg) 

-- 

Iberall 

(analytic) 

-- 

-- 

9.940 

1.381 

2.095 

-- 

; Watts 
(exper * ) 


-- 

1 

1 

-- 

-- 

i 

UPD: 

continuity 

-- 

-- 

9.942 

1.232 

10.52 

2 

UPD: 

momentum 

-- 

i 

i 

9.942 

1.232 

10.52 

3 

Equilibrium: 

continuity 

Increments 

per 

cycle 

50 

9.942 

1.332 

2.81 

100 

9.943 

1.336 

1.68 

4 

Equilibrium: 

momentum 

Increments 

per 

cycle 

50 

9.942 

1.332 

2.81 

100 

9.943 

1.336 

1.67 

5 

Implicit: 
no pressure 
linking 

Increments 

per 

cycle 

200 

9.943 

1.338 

1.10 

6 

UPD: 

n - step 
implicit 

Increments 

per 

time step 

2 

9.942 

1.302 

6.06 

3 

9.942 

1.318 

3.78 

5 

9.942 

1.329 

2.86 

10 

9.942 

1.335 

1.72 

7 

n - step 
globally 
implicit 

Increments 

per 

cycle 

5 

9.942 

1.052 

13.73 

10 

9.942 

1.228 

10.66 

20 

9.942 

1.304 

6.15 

8 

n * step 
implicit : 
character- 
istically 
determined 
time step 

Increments 

per 

time step 

2 

9.938 

1.301 

6.18 

3 

9.938 

1.318 

3.74 

5 

9.938 

1.329 

2.85 

10 

9.938 

1.335 

1.69 
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Table A. 4 6 inch transmission line with 400 Hz excitation 


Reference 

Number 

Analysis 

Temporal 

Discretisation 

Parameter 

Parameter 

Value 

Char. 

Number 

Amplitude 

Ratio 

Phase 

Angle 

(deg) 


Iberall 

(analytic) 

i 

i 

-- 

5.716 

3.885 

10.886 

-- 

Watts 
(exper . ) 

-- 

-- 

l 

1 

5.6 

-- 

i 

UPD: 

continuity 


-- 

5.718 

1.208 

40.12 

2 

UPD: 

momentum 

-- 

-- 

5,718 

1.209 

40.18 

3 

Equilibrium: 

continuity 

Increments 

per 

cycle 

50 

5.721 

3.180 

15.86 

100 

5.721 

3.301 

9.04 

4 

Equilibrium: 

momentum 

Increments 

per 

cycle 

50 

5.720 

3.180 

15.85 

100 

5.721 

3.301 

9.03 

5 

Implicit : 
no pressure 
linking 

Increments 

per 

cycle 

100 

5.721 

3.301 

9.07 

6 

UPD: 

n - step 
implicit 

Increments 

per 

time step 

2 

5.719 

1.997 

39.21 

3 

5.720 

2.470 

32.55 

9 

5.721 

3.171 

14.76 

17 

5.721 

3.285 

8.96 

7 

n - step 
globally 
implicit 

Increments 

per 

cycle 

5 

5.717 

1.029 

38.78 

10 

5.718 

1.791 

42.70 

20 

5.720 

2.604 

31.52 

8 

rx - step 
implicit : 
character- 
istically 
determined 
time step 

Increments 

per 

time step 

2 

5.712 

1.972 

i 

38.84 

3 

5.712 

2.412 

34.08 

9 

5.711 

3.123 

19.48 

17 

5.711 

3.281 

14.80 






















































































Table A. 5 24 inch transmission line with 55 Hz excitation 


Reference 

Number 

Analysis 

Temporal 

Discretisation 

Parameter 

Parameter 

Value 

Char. 

Number 

Amplitude 

Ratio 

Phase 

Angle 

(deg) 

-- 

Iberall 

(analytic) 

-- 

-- 

10.340 

1.299 

3.607 

-- 

Watts 
(exper . ) 

-- 

-- 


-- 

-- 

i ; 

UPD: 

continuity 

-- 

-- 

10.341 

1.162 

8.62 

2 ! 

UPD: 

momentum 


-- 

10 . 341 

1.162 

8.62 

3 

Equilibrium: 

continuity 

Increments 

per 

cycle 

50 

10.341 

1.234 

3.22 

100 

10.341 

1.237 

2.50 

4 

Equilibrium: 

momentum 

Increments 

per 

cycle 

50 

10.341 

1.234 

3.23 

100 

10.341 

1.237 

2.50 

5 

Implicit: 
no pressure 
linking 

Increments 

per 

cycle 

170 

10.341 

1.238 

2.23 

6 

UPD: 

n - step 
implicit 

Increments 

per 

time step 

2 

10.341 

1.214 

5.47 

3 

10.341 

1.226 

4.03 

5 

10.341 

1.233 

3.32 

10 

10.341 

1.237 

2.53 

7 

n - step 
globally 
implicit 

Increments 

per 

cycle 

5 

10.341 

1.032 

11.30 

10 

10.341 

1.161 

8.68 

20 

10.341 

1.214 

5.50 

8 

n - step 
implicit : 
character- 
istically 
determined 
time step 

Increments 

per 

time step 

2 

10.338 

1.215 

5.16 

3 

10.338 

1.226 

4.25 

5 

10.338 

1.234 

3.50 

10 

10.338 

1.238 

2.84 
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Table A. 6 24 inch transmission line with 110 Hz excitation 


Reference 

Number 

Analysis 

Temporal 

Discretisation 

Parameter 

Parameter 

Value 

Char. 

Number 

Amplitude 

Ratio 

Phase 

Angle 

(deg) 

-- 

Iberall 

(analytic) 

-- 

-- 

5.170 

4.252 

25.419 

-- 

Watts 
(exper . ) 

-- 

-- 

-- 

4.8 

-- 

1 

UPD: 

continuity 

— 

-- 

5.171 

1.008 

40.10 

2 

UPD: 

momentum 

-- 

-- 

5.171 

1.008 

40.11 

3 

Equilibrium: 

continuity 

Increments 

per 

cycle 

50 

5.174 

3.070 

19.60 

100 

5.174 

3.209 

13.26 

4 

Equilibrium: 

momentum 

Increments 

per 

cycle 

50 

5.174 

3.070 

19.60 

100 

5.174 

3.209 

13.26 

5 

Implicit : 
no pressure 
linking 

Increments 

per 

cycle 

90 

5.174 

3.196 

13.99 

6 

UPD: 

n - step 
implicit 

Increments 

per 

time step 

2 

5.172 

1.743 

43.62 

3 

5.173 

2.211 

38.82 

10 

5.174 

3.063 

19.55 

20 

5.174 

3.204 

13.30 

7 

n - step 
globally 
implicit 

Increments 

per 

cycle 

5 

5.171 

1.007 

39.91 

10 

5.172 

1.743 

43.97 

20 

5.173 

2.512 

33.88 

8 

n - step 
implicit: 
character- 
istically 
determined 
time step 

Increments 

per 

time step 

2 

5.167 

1.798 

42.56 

3 

5.167 

2.262 

37.81 

10 

5.166 

3.069 

19.14 

20 

5.166 

2.942 

16.79 
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A. 7 CONCLUSION 


The results confirm that the discretised primitive integral conservation 
balances have the following features: 

they are mathematically symmetrical with respect to continuity and 
momentum equation linking. 

when they are solved implicitly, as the time increment is decreased, the 
pressure field solution obtained using either a continuity or momentum 
linked equation set approaches the solution obtained explicitly from the 
continuity equation via an equation of state. 

in terms of the transmission line, their solution is independent of the 
form of the implicit numerical integration algorithm used. 

Therefore, in the light of these confirmations, it seems reasonable to 
conclude that the discretised pressure or mass flux linked primitive 
conservation balances used In the simulations are not unconditionally valid 
for modelling transmission line information propagation effects at 
characteristic numbers less than about 24. To the extent that these specific 
effects are also present in Stirling machines (in particular, the SPDE and 
SPRE) , it is consistent to infer t hat Stir ling machine simulations based on 
these primitive conservation balances are also open to question. 

In particular, at low characteristic numbers, neither the pressure - 
linked nor mass-linked formulations appear to replicate the physical behavior 
described by the Kirchoff momentum equation or by the wave equation, despite 
attempts to impose a characteristics solution on the pressure-linked 
formulation using either a UPD or a characteristically determined, variable 
time step hypothesis. 

This position is supported by Iberall's analysis itself, by the results 
of Organ's linear wave equation approach (OJ89) as well as by MacCormack 's 
analysis (Mc82), . 
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PART £ 


ENTROPY TRANSPORT 

i 


B.l INTRODUCTION 

The Investigation of entropy transport resolved Itself into two tasks. 

In the first task, a previously developed suite of codes used for simulating 
the Mechanical Engineering Test Rig (METR) (Go90) was applied to Kurzweg's 
apparatus and validated against his experimental and closed- form analytic 
data. Thereafter, an entropy transport equation was developed and evaluated 
in the validated codes. These activities may be resolved into the following 
specific tasks: 

a. The adaptation of the Mechanical Engineering Test Rig (METR) one- and 
two-dimensional codes (Go90) to Kurzweg's enhanced diffusion oscillating 
flow test apparatus (KZ84) . 

b. The formulation of an entropy postulate and its rigorous transformation 
into an entropy transport equation to be tested in a code applied to 
Kurzweg ' s apparatus . 


B.2 INITIAL OBSERVATIONS 

Kurzweg's apparatus is based on incompressible (water) fluid flow. 
However, in view of NASA's interest in the compressible oscillating flows 
occurring in Stirling cycle machines, it was originally intended to simulate 
Kurzweg's apparatus using a compressible fluid. This was also logistically 
prudent since the METR codes are strictly appropriate for compressible flows. 
However, early on during task a it became clear that the original Intention of 
using a compressible fluid was not feasible for at least the following 
reasons : 

reference KZ84 provides insufficient geometrical detail about the 
variable volume cavities at either end of the capillary tube bundle to 
enable the compressible flow boundary conditions to be specified 
completely . 

in order to remain within the laminar flow regime specified for the 
applicability of the analysis described in KZ84, very small compressible 
flow oscillation frequencies and tidal displacements are required. Under 
these conditions, other effects not accounted for by Kurzweg (such as 
buoyancy effects) may be expected to influence the results. In this 
case, the validity of the experimental data in KZ84 under compressible 
flow conditions is questionable. 
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Thus it was decided to enhance the METR codes to describe incompressible 
flows also so that they could be applied directly to Kurweg's apparatus 
without the above two limitations. 

Initial validation efforts with the two-dimensional code at the chosen 
test point (see section B.3), consistently yielded effective dif fusivities 
four to five times larger than those ^measured or predicted analytically. The 
simulation was found to be very sens itive t o axial discretisation so that each 
decrease in discrete volume axial length decreased the simulated effective 
diffusivity, edging it ever closer to the nominal experimental value. At a 
discretisation of 1 discrete volume per millimeter, the overprediction had 
been reduced to a factor of 2 at the e xpense of a very large increase in 
computational effort, almost overwhelming the computer equipment available for 
the project. 

These observations are prima facie evidence of the ubiquitous "false 
diffusion" problem that plagues numerical analysis in general and integral 
analysis in particular. Hence further progress mandated the development of an 
algorithm for eliminating the false diffusion problem inherent in the 
discretised integral description of the continuum mechanics used in the 
simulation. This made inclusion of entropy transport in the simulations of 
vital importance, since false diffusion correction schemes do have the 
potential for violating the second law of thermodynamics. 

In order to develop such a false diffusion correction algorithm, it was 
judged prudent to do the development work using a one -dimensional simulation 
of Kurzweg's apparatus rather than a two-dimensional code since this was felt 
to be more time -effective in the long run. Furthermore, since the one- 
dimensional simulation uses the standard Kays and London (KL64) steady- state 
heat transfer coefficient and friction factor correlations, such a simulation 
provides another experimentally based test point for evaluating the accuracy 
of such correlations in describing oscillating flows. 


B.3 TEST CONDITIONS 

The test apparatus simulated is shown in figure B.l. Essentially, the 
device consists of an upper and a lower cavity joined by an acrylic tube 

containing a bundle of glass tubes. These tubes apparently are not fitted 

Into a manifold at either end so that the net flow area consists of the 
internal tube areas as well as the triangular areas between the tubes (figure 
B.2). The external axial faces of the hot and cold reservoirs are fitted with 
flexible membranes so that when the entire apparatus is oscillated by a 
shaker, the fluid (water) will undergo corresponding oscillations. Finally, 
the apparatus is cocooned in insulation so maintaining a nominally adiabatic 
boundary with the surroundings. 
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Figure B.l Kurzweg's test apparatus (reproduced from KZ84) 



triangular flow 
area 


Figure B.2 Typical segment of the tube bundle cross-section 


The experimental methodology employed is open to criticism. Initially, 
the lower (cold) cavity is filled with dyed cold water, while the tube bundle 
and upper cavity are filled with hot water. Hence, the tidal displacement is 
determined by measuring the extent of penetration of the dyed water into the 
tube bundle. It must be noted that such simplicity is prone to errors such 
as ; 


the mass and thermal diffusion of the cold water into the hot water 
which leads to a "fuzzy" boundary zone making the precise measurement of 
tidal displacement difficult. 

radial non-uniformity of the flow oscillations arising from the specific 
geometry of the shaker/flexible membrane interface, that is, the tidal 
displacements at the center of the tube bundle may be different from 
those at the containing wall. 
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surface tension effects within the capillaries which would complicate 
the external measurement of the mid-meniscus position. 

Kurzweg and Zhao do not appear to acknowledge these problems in the 
light of the good agreement they obtain (in logarithmic terms) between their 
analytic and experimental results (KZ84, figure 2). However, since they use a 
data fitting approach to obtaining this agreement (by adjusting the Prandtl 
number via the choice of system temperature) , it is reasonable to suppose that 
many of the systematic errors which would degrade this good agreement are 
buried in the data fitting process. 

The upper and lower cavity temperatures are measured via thermometers 
read at one minute intervals for a total experimental run time of 6 minutes. 
The heat flowing through the tube and resultant effective diffusivity are 
determined exclusively from these measurements. In particular, the resulting 
inferred temporal temperature gradient of the cold cavity (from which the net 
heat flowing down the tube is computed) is highly approximate in corrifarison 
with the instantaneous axial temperature gradient used in the effective 
diffusivity calculation (KZ84, equation 1) (a time-averaged axial temperature 
gradient would have been more appropriate) . Hence this additional source of 
systematic error is also embedded in the Prandtl number data fit. 

Needless to say, these issues complicate the simulation, particularly 
when the objective is to obtain agreement with the experimental and 
"validated" analytic results. However, of perhaps even more importance from a 
simulation perspective, is the difficulty of simulating the initial 
discontinuous (literally "brick wall") temperature profile at the tube bundle 
/ cold cavity interface. Such a condition has been documented in the 
literature (for example, BB75) to magnify the effects of false diffusion in 
producing erroneous advective fluxes. 

Hence within these limitations, the following simulation approach is 
adopted . 

a. All the available physical parameters describing the apparatus and test 
conditions given in KZ84 are used as described in table B.l. 

b. The three additional geometrical parameters needed to implement the 
simulation were determined by direct measurement of figure B.l in the 
hope that the relative scale of the figure is somewhat representative of 
reality. These parameters are: 

Hot/cold cavity cross -sectional area ratio = 1 
Stationary hot/cold cavity length ratio = 1.714 
Cold cavity / casing cross-sectional area ratio =-4.8 

Fortunately, under incompressible flow conditions, errors in these 
ratios do not effect the flow within the tube bundle provided mass flux 
boundary conditions are used. They do have an effect when pressure 
boundary conditions are invoked. 
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Table B.l Explicit apparatus and test parameters 


Parameter Description 

Parameter Value 

Number of capillary tubes 

31 

Capillary tube length 

200 mm 

Tube bundle casing diameter 

12.7 mm 

Capillary tube internal diameter 

1 mm 

Net flow area 

0.67 cm 2 

Cold cavity volume 

114 cm 3 

Cold cavity initial temperature 

22 "C 

Hot cavity initial temperature 

78 °C 


c . 


The test point was chosen from KZ84 figure 2 (a logarithmic plot) in 
terms of Kurzweg's correlation factor %% which is defined as: 


K % * Ax 


2 yw 


(B.3.1) 


where : 


Ax = tidal displacement 
v = kinematic viscosity 
w - angular velocity 
r - capillary tube radius 


The nominal chosen value of k K **= 122 cm z / s corresponds to the 
middle of the test range where evidently three tests were carried out at 
the same K K (or three data points from the same test) . These three 
points yield values just above that calculated using the analytic result 
given by equation (15) in KZ84 , Furthermore, at this point in the 
correlation, the analytic result is bounded below by two additional 
experimental values at a slightly greater (in a logarithmic context) k k 
of 149 cm 2 /s . Hence the "fit” of the analytic result appears to be very 
good at this point. The test K K is achieved by choosing an oscillation 
frequency of 6 Hz and a tidal displacement of 37.83 mm which both fall 
within the parameter ranges actually tested. 


d. As noted by Kurzweg and Zhao, the accuracy of their analysis is 

critically dependent upon the choice of temperature at which the fluid 
parameters are evaluated, particularly the Prandtl number, which is 
strongly dependent on temperature in the 22 to 78 °C range. Therefore, 
the fluid parameters used in the simulations are chosen to be those at 
60 °C as recommended in KZ84 (figure 2). The relevant parameters for 
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water at this temperature together with those of glass (necessary for 
modelling the capillary tube walls) are given in table B.2. 


Table B.2 Fluid and wall material properties 


Property 

Water (Sc79) 

Glass (ASH81) 

Density (kg/m 3 ) 

983.3 

2470 

Specific heat capacity (J/kg.K) 

4191 

750 

Thermal conductivity (W/m.K) 

0.65 

1.0 

Dynamic viscosity (kg/m.s) 

47.0 x 10' 5 

— 

Prandtl number 

3.01 



It must be pointed out that a better approach would be to include these 
parameters in the simulation as temperature dependent entities. 

However, in the interest of reducing the differences between the 
simulation and the analysis (the major objective being to use the 
analysis as fitted to the experimental data as the validation standard), 
it was decided to proceed on the above basis. 


e . 


The method used to calculate the effective diffusivity in the simulation 
corresponds to that of KZ84 and is given by: 


. fll Vc ^ c 

[ A jtuJbe 


(B.3.2) 


In applying equation (B.3.2), the rate of temperature change In the cold 
cavity is computed as the net cyclic temperature change rate, that is, 
the difference between the cold cavity temperatures at the end and 
beginning of a cycle divided by the period. This enables the change In 
effective diffusivity to be tracked as a function of time so that the 
simulation can be' terminated as soon as the rate of change of the 
simulated effective diffusivity becomes small. Ultimately, of course, 
the effective diffusivity becomes infinite when the hot and cold 
cavities reach thermal equilibrium. This achieves the same physical 
purpose as Kurzweg's experimental procedure while minimizing the amount 
of computation. In practice, a simulation period of .5 minutes proved 
more than adequate under these conditions. 


B .4 THE FALSE DIFFUSION CORRECTION METHODOLOGY 

False diffusion arises in transient continuum mechanics numerical 
analysis in the discretisation of the advection terms in the conservation 
balances. Typically, when using any variation of upwind difference or upwind 
parameter methodology (whether the second upwind difference formulation of 
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Gentry, Martin and Daly (GM66) or the equilibrium analytic approach of 
Patankar (Pa80)), the advected flux on the control volume boundary is biased 
towards the upwind advected parameter rather than the actual value of the 
parameter at the boundary. In general, the actual value of the advected 
parameter Is different from both the upwind and downwind values In a 
continuous system. In Stirling machine analysis, false diffusion is known to 
be a significant cause of error in modelling the working fluid in the 
regenerator (originally pointed out by Gedeon (Ge84)). 

The literature on false diffusion is quite voluminous as it long has 
been recognized as one of the principal difficulties in numerical analysis. 
However, most of the algorithms developed for correcting or eliminating false 
diffusion are based on discretisations of the differential conservation 
balances, particular examples being the SHASTA algorithm of Book, Boris and 
Hain (BB75) and the sequential backward difference predictor / forward 
difference corrector methodology of MacCormack (Mc82). Unfortunately, none of 
these schemes works well in an integral environment where there are no 
parameter gradients to be exploited in the advection terms. 

The first attempt made was to invoke the upwind linear extrapolation 
method successfully used in Go87 for modelling false diffusion in the 
regenerator of Stirling machines. Theoretically, such an approach could be 
expected to work well in simulating Kurzweg's apparatus, since in time, the 
temperature profile between the hot and cold cavities becomes linear. 

However, at startup, in the presence of the discontinuous temperature gradient 
at the cold cavity entrance, the method fails badly, and later, even under 
milder temperature gradient conditions, the method yields an over-prediction 
of the effective diffusivity. 

Hence, after much analysis and computation, an apparently effective and 
unconditionally stable false diffusion correction methodology was developed. 
The methodology may be described in terras of the following generalized 
algorithmic sequence: 

a. Predict the parameter field using a conventional second upwind 

difference spatial discretisation (Go87). In general, an explicit, 

implicit or hybrid integration scheme may be used for this step. 

b. Calculate the source term from: 
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(B.4.1) 


K source 




where : 

fuw = tl -1 if 9ni * ° 

^ uw = If £7ni ^ 0 

V£ « f(^i) 

It may be noted that no restrictions need be placed on the form of the 
interpolating function f provided it is bounded symmetrically , for 
example, f may be a polynomial of odd order. 

The correction field A^> is found implicitly from: 

= -ii)dA + K source (B.4.2 


where : 

Wuw - A V>i- 1 if ~9ni * 0 
M>uw s A V>i if -9n± < 0 


d . Calculate the corrected parameter field ^ given by: 


(B.4.3) 


e. Iterate from step a until the change in the corrected parameter field V 

becomes appropriately small. 

In practice, the iteration required in step e does not add to the 
iterations required in any case to implement the pressure or non-pressure 
linked algorithms discussed in part A. Two interpolating functions have been 
tested for use in step b, namely linear and cubic polynomials. As shown by 
the results (section B.7), both schemes produce almost the same results, 
however, during the transient start-up phase, the cubic function yields a 
better approximation of the temperature field at the cold cavity entrance. 

A negative consequence of the f a 1 s_e_ dif fusion correction algorithm 
described above is that it necessitates one more matrix inversion per 
conservation equation per iteration. However, in terms of an integral 
analysis, this is more than offset by at least a six-fold reduction in the 
required discretisation combined with better accuracy at the coarser 
discretisation than obtained using the upwind approach alone at the finer 
discretisation. In this context, it must be pointed out that the flux 
correction methodology developed is still embryonic and is capable of much 
enhancement . 
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B . 5 DEVELOPMENT OF AN ENTROPY TRANSPORT EQUATION 


The approach adopted in deriving the entropy transport equation is 
identical in nature to that documented in Go87 for the mass, momentum and 
energy balances. The essential difficulty (at least philosophically) is in 
defining an appropriate entropy postulate. One traditional approach is to 
define entropy as being that which requires thermodynamic processes to be 
unidirectional, for example, according to Wark (Wa77): 

"Any system having certain specified constraints and having an 
upper bound in volume can reach from any initial state a stable 
equilibrium state with no net effect on the environment." 

The corollary to this postulate, known as the Kelvin-Planck statement (Wa77) , 
is sufficient to allow proof of Clausius' inequality as a theorem, that is for 
a closed system: 

f*g:SO (B.5.1) 

J T 


which in turn allows the "increase in entropy" principle to be derived, or. 


dH> I£ 

T 


(B.5.2) 


where H denotes entropy. 

Equation (B.5.2) is invoked by some authors directly as an entropic 
postulational basis, for instance Slattery's entropy postulate (S181) that: 

"The minimum rate of production of entropy in a body is 
proportional to the rate of energy transmission to the body." 

or symbolically: 


d 

Si 


[v lm f~ AV> Ja^/t * 


-n)dA 



fldv 

T 


(B.5.3) 


where the terms on the RHS represent the entropy changes produced by contact 
energy transmission to the body through its bounding surface and by external 
and mutual energy transmission respectively. 

However, from a computational point of view, neither of these postulates 
is determinate, that is, they do not allow entropy to be definitely 
calculated. Hence for this reason, the following entropy postulate based on a 
statement originally proposed by Truesdell and Toupin (TT60) as the basic 
assumption of thermodynamics is preferred: 
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Postulate V A dimensionally independent scalar parameter S (termed specific 

entropy) and the substate densities are sufficient to determine 
the specific internal energy of a material particle independently 
of time, place, motion and stress. 


That is : 

d-ofi. P± Bl.x- p) 


(B.5.4) 


This postulate is not equivalent to those quoted above in the sense that 
it is not based upon a manifestation of entropy, but rather is an equation of 
state for internal energy which requires that entropy be a fundamental 
property. However, the consequences of equation (B.5.4) must satisfy the 
previous conventional postulates as expressed in a general fashion by equation 
(B.5.3). 


When a particular body is thermodynamically homogeneous, that is every 
particle p is behaviorally identical, equation (B.5.4) may be simplified to: 


6 = 6 



V, 



(B.5.5) 


where T is the mass (or mole) fraction of species i. In terms of the single 
component working fluids being considered here (and hence also for Stirling 
machines in most cases) , it is useful for the sake of simplicity to proceed by 
considering single specie entities only, that is: 


V = u{ s, v) 


(B.5.6) 


Differentiating and noting that 



and Ps 



yields : 


dU = Tdk - PdV 


(B.5.7) 


Taking the substantive derivative, multiplying through by p, substituting Go87 
equation (C. 7) (which relates the change in size of a discrete volume to the 
motion of its enclosing boundaries) and simplifying: 


p ^ -p(v*v) 


T5t 


(B.5.8) 


The substantive differential energy balance at every point in a material body 
is given by (Go87) : 


p^£ = pi + (r:Vv) -P(V*v) -V.$ 


(B.5.9) 
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Equating (B.5.8) and (B.5.9), expanding and rearranging: 


D5 , p£ + (r:Vv) i vr (B.5.10) 

T5t T T T 


Integrating (B.5.10) over the material body volume V im and substituting Go87 
equations (C.9) and (C.10) yields after simplification and rearrangement: 


d 

lit 


f flSdV = f 

fpE + (T: Vv) 

- JiVT Idv + f 

? 

— » 

(m)\ T T 

r 2 J 1 


(B.5.11) 


A («) 


Rearranging the right hand side: 


f pE dV = 1 

( 4 ; • -n)dA + | 

£E dv + f 

f(r:Vv) 

1 

k 

< 

(B.5.12) 

jV(m) J 

Am) T J 

V) T >< 

m)( T 

X 2 J 



Comparing equations (B.5.12) and (B.5.3) and noting that the last 

integral on the RHS of (B.5.12) is strictly positive (because £= -kVT ), 

reveals that the postulational basis of equation (B.5.4) indeed satisfies the 
conventional entropy postulate as defined by equation (B.5.3) while 
simultaneously being determinate. 


Converting equation (B.5.10) into partial temporal derivative form and 
substituting into the generalized transport equation (Go87, equation (C . 12) ) 
yields the final version of the entropy balance: 


d 


m EM 

At 




pE + (T: Vv) 




+ 



v (S) ).-n}dA 


(B.5.13) 


In Implementing the entropy transport equation, it was found that great 
care must be taken in locating the entropy boundary (within which the cyclic 
entropy is computed) and in describing the thermal transport across that 
boundary. In this context, it should be noted that, by definition, the second 
law of thermodynamics treats isolated systems so that Clausius' inequality 
applies over the integral of the subsystems of which the isolated system is 
comprised. In the limit, the only true isolated system is the entire 
universe, so that extracting a computationally isolated system which can in 
general exchange heat with its superordinate body is not obvious in all cases. 

The approach adopted may be visualized in terms of figure B.3. 
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Figure B.3 Isolating boundary entropy contribution 

Consider a segment of the isolating boundary (& n )± 3 across which flows a mass 
flux (g n )± 3 and a thermal diffusion flux ($ n )i s . The operating entropic 
principle governing the isolating boundary is: 

The net contribution of the entropic and thermal fluxes crossing the 
isolating boundary is assigned to the isolated system entropy 
irrespective of the spatial location of these fluxes (which in general 
span a part of the isolated system as well as the bounding environment) * 

Hence, applying equation (B.5.13) to the situation depicted in figure B.3 
yields : 


d m Etf _ 

d m EM 

+ r : 

{(4 r n)is 7r Jb} Hv 

dt 

dt 

. j 

internal 

T % 


+ r 

-[(Qn)i3^ T sys) HU 

* f 

(^)is 



' 2 

T ays 

JtA.Ji, 

r sys ~ T b 


(B.5.14) 


where a is an appropriate fraction for delineating the integration sub-volume 
(typically a - .5 for a conformal Cartesian mesh). It should be noted that 
the advection entropy flux does not appear in equation (B.5.14) because its 
net contribution to the isolated system entropy is zero. 

The success in using the full entropy transport equation including 
generalized Isolated boundary heat transport may be judged from the results 
presented in section B.7. As envisaged when proposing this study, Kurzweg's 
laminar flow apparatus provides an ideal test bed for such an entropy 
transport equation, since the transport of entropy is clearly visible from a 
system perspective and not compartmentalized between system components as 
occurs, for example, in Stirling machines, 

B. 6 OTHER ISSUES 

During the process of testing the false diffusion correction algorithm, 
it was found that careful attention had to be paid to modelling the conduction 
heat transfer in the capillary tube walls. In particular, the effective 
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diffusivity was found to be sensitive to the manner in which the wall/fluid 
interface is discretised, both in the capillary tubes as well as in the hot 
and cold cavities. Although in the simulations (both one- and two- 
dimensional) the capillary tube walls are modelled as axially discretised, 
one -dimensional entities and this approach seems to yield good results, there 
is a suspicion that a one -dimensional approach is probably inadequate, 
particularly at high radial heat transfer rates. Thus, in the future, thought 
should be given to modelling Stirling heat exchanger walls as radially 
discretised two-dimensional entities to determine whether the above suspicion 
has any basis. 

Finally, as the results show, the limitations of a first order temporal 
integration scheme became manifest, particularly in the two-dimensional 
simulation. However, owing to time constraints, the first order temporal 
integration algorithm in the METR two-dimensional code could not be upgraded 
to second order (a second order one -dimensional code had been developed and 
tested previously (Go90)), 


B . 7 RESULTS 

The results of all the simulations undertaken are summarized in table 
B.3. At the chosen test point, the effective diffusivity produced by 
Kurzweg's analysis amounts to 4.8 cm z /s with the corresponding experimental 
values ranging between 5.8 and 7.4 cm z /s . Using an axial discretisation of 30 
discrete volumes (.15 discrete volumes per millimeter) and a temporal 
discretisation of 100 increments per cycle, the first order one -dimensional 
simulation with no false diffusion correction over-predicted the effective 
diffusivity by a range of 2.9 to 4.4 times (corresponding to the upper 
experimental and analytic values respectively). As noted above, including the 
false diffusion correction produced a simulated value for effective 
diffusivity within the experimental range and at most 36% greater than that 
analytically predicted. In contrast, upgrading the one -dimensional simulation 
to second order temporal accuracy, reduced the over-prediction of effective 
diffusivity with no correction to a range of 2.6 to 4.1 times. Including 
either linear or cubic false diffusion correction produced a result within 
3.5% of the analytic value. 

Applying the second-order simulation at two other values of K% using 
cubic false diffusion correction yields the results compiled in table B.4. 

As discussed in section B.3 above, in view of the experimental 
uncertainties, perhaps the "fitted" analytic value Is Indeed a reasonable 
estimate of the actual effective diffusivity, in which case, discrepancies 
less than 3,5% are an adequate validation of the false diffusion correction 
methodology , 
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Table B.3 Simulation results 
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Table B.4 Effective diffusivity comparison 


Correlation 
Factor Xjj 
( cm 2 /s) 

Effective Diffusivity (cm z /s) 

Analysis/ 

Simulation 

Kurzweg' s 
Analysis 

Experimental 

1-d Simulation 
(cubic correction) 

Discrepancy 

(%) 

122.0 

4.8 

5.8 - 7.4 

4.9459 

3.04 

286.1 

11.31 

11.1 

11.0539 

2.26 

818.5 

32.36 

19.1 - 25.9 

31.2692 

3.37 


However, of perhaps more interest, is the observation that all the one- 
dimensional simulations were carried out using the nominal steady-state Kays 
and London correlations without the necessity of any heat transfer coefficient 
multipliers. For all these simulations, the Valensi number (square of the 
Womersley number) is 19.7 which remains constant because it is independent of 
tidal displacement. was varied in the simulations via the tidal 

displacement while the frequency was kept constant at 6 Hz. 

At a Valensi number of about 20, Uchida's analysis (Uc56) of oscillating 
laminar flows as quoted by Simon and Seume (SS86), yields amplitude 
coefficients of pressure drop and wall shear stress of about 3.1 and 1.3 
respectively. Also, at this Valensi number, Uchida's analysis predicts lead 
phase angles of pressure drop and wall shear stress of approximately 68° and 
28° respectively. Noting that for steady flow the amplitude coefficients are 
unity and the lead phase angles are zero, it is evident that the simulations 
represent a flow regime which is different from steady flow (particularly in 
phase lead terms). Hence in view of the good agreement obtained in table B.4, 
an argument can be made for the adequacy of the steady- state correlations for 
predicting heat transfer under laminar oscillating flow conditions at Valensi 
numbers less than 20 when the effects of false diffusion are minimized (if not 
eliminated) . 

In developing the two-dimensional results, the importance of temporal 
discretisation has already been alluded to. For this, reason, the high 
temporal discretisations required to minimize the limitations of the first 
order temporal integration algorithm used placed a severe load on the 
available computing resources which constrained the number of runs carried out 
to the half -minute closure point. 

The axial discretisation was maintained at 30 discrete volumes while the 
radial discretisation could be comfortably limited to 8 discrete volumes. 

With the false diffusion correction activated, increasing either of these 
values produced a negligible change in the end result. 
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Hence in this context, other than the temporal discretisation, the only 
significant parameter was found to be the mass flux boundary condition. Two 
conditions were tested, namely; 

specifying a radially uniform flow at the cold end tube entrance equal 
to the rate of change of the cold cavity volume. 

using the rate of cold cavity volume change to generate a sinusoidally 
varying pressure which in turn allows the entrance mass fluxes to be 
calculated (a somewhat nebulous approach in the light of the paucity of 
data describing the hot and cold cavities). 

In both cases, the hot cavity is treated as being isobaric at 
atmospheric pressure and gravity induced buoyancy forces are included in the 
momentum balance . 

In all cases, the uniform flux boundary condition produced slightly 
better results as might be expected. A shown in table B.3, only the linear 
false diffusion correction methodology was tested. At a temporal 
discretisation of 400 increments per cycle, the uniform flow boundary 
condition simulation produced an effective diffusivity between the 
experimental and analytic values, while increasing the discretisation to 600 
increments per cycle yields a simulated diffusivity 6.5% higher than the 
analytic result and shows the diffusivity converging to the value predicted 
using the second order one -dimensional code. Nevertheless, while yielding a 
satisfactory result, a discretisation of 600 increments per cycle is not 
practical for an implicit code. 

At the half minute cut-off, the false diffusion correction was 
deactivated with a resultant escalation of the effective diffusivity to an 
unconverged value in excess of 30 cm 2 /s. This is indicative qualitatively of 
the effectiveness of the false diffusion correction algorithm in two 
dimensions. 

The results produced by the entropy transport equation are depicted in 
figures B.4 through B.9. Figure B.4 shows that the entropy is highest midway 
along the tube bundle and decreases towards both cavities. The hot cavity 
always has a lower entropy than the cold cavity. Numerical output from the 
simulation (not shown) confirms that the cyclic integral of the entropy over 
the entire fluid system increases monotonically . 

The cyclic temperature surface is depicted in figure B.5. Essentially, 
this shows a fairly invariant linear temperature profile flattening 
alternately in the regions adjacent to the hot and cold cavities as the mass 
flux reverses direction. 

Combining figures B.4 and B.5 yields the temperature - specif ic entropy 
(T-H) diagrams of figures B.6 to B.9. Within the tube, figures B.7 and B.8 
show closed cycles whose areas are exactly equal to the generated cyclic 
irreversibilities at those locations. As the tube acts as a thermal transfer 
device between the hot and cold cavities and is externally adiabatic, closure 
of the T-S diagrams is thermodynamically necessary and thus an indication that 
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the entropy transport equation is operating correctly. All the entropy 
generated in the tube is pumped towards the hot and cold cavities, both of 
whose entropies must increase monotonically which indeed is the case as shown 
in figures B.6 and B.9. In figure B.6, during the first part of the cycle 
when the flow is positive (that is, from the cold to the hot cavity), the 
specific entropy increases marginally while the temperature drops (dissipation 
Is not explicitly included in the cavities because of the absence of adequate 
geometrical data, see section B.3). When the flow reverses, the entropy 
generated In the capillary tubes is transported into the cold cavity so 
producing the requisite monotonic cyclic entropy increase. 

The inverse process occurs In the hot cavity as shown in figure B.9, 
During the first half of the cycle when the flow is positive, entropy is 
transported into the hot space while Its temperature decreases. Notice that 
over the cycle, the hot space has cooled while the cold space has become 
warmer. 

At first sight, It may seem unreasonable that the temperature of the 
cold cavity decreases during positive flow and the hot cavity temperature 
increases during negative flow. However, bearing In mind that in an integral 
analysis the reported temperatures are the volume -averaged temperatures and 
that an axial temperature gradient exists through the cavities, then as mass 
is removed from a cavity, its volume - averaged temperature must change. Hence, 
in the cold cavity, the removed fluid is hotter than the remaining fluid and 
hence the volume -averaged temperature decreases, with the inverse process 
occurring In the hot cavity. This Is a direct consequence of the flux 
correction methodology and reveals how it provides a more accurate description 
of the continuum mechanics than the upwind difference approach. In 
particular, in the absence of diffusion heat transfer from the ends of the 
capillary tubes, the upwind algorithm yields constant temperatures in the 
cavities during their respective periods of mass efflux. 
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Figure B. 4 ENTROPY SURFACE 

Kurzweg’s Oscillating Flow Apparatus 

1 —Dimensional Simulation 



39 


Note: Distance measured from cold end 








Figure B.5 TEMPERATURE SURFACE 
Kurzweg’s Oscillating Flow Apparatus 
1 —Dimensional Simulation 
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Note: Distance measured from cold end 






Figure B.7 TUBE 1 
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B.8 QO ^ClVSm 

The significant results obtained from the work carried out in part B, 
have been the development of a tentatively effective false diffusion 
correction methodology, its demonstration in one- and two-dimensional 
simulations and the development and evaluation of a rigorous entropy transport 
equation. Using Kurzweg's oscillating flow apparatus as the target hardware, 
in terms of the prevailing geometrical and experimental uncertainties, the 
simulated data compare favorably with those generated analytically and 
experimentally . 

A comparison of the simulation and experimental/analytic data suggests 
that under laminar, incompressible flow conditions at Valensi numbers less 
than 20, standard steady- state heat transfer correlations are adequate for 
predicting oscillating flow boundary heat transfers provided that the effects 
of false diffusion are accounted for in the simulation. 

Use of a full entropy transport equation in the simulation provides 
useful physical insights into the irreversibility processes occurring within 
oscillating flows. 

The results produced should shed some additional light on the loss 
mechanisms occurring within Stirling cycle machines, at least in terms of 
isolating possible sources of error in existing simulation design codes. 


j 
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CLOSURE 


The work carried out in part A suggests that the next logical steps for 
investigating information propagation effects in Stirling machines might be: 

the development and validation of a simulation based on a solution of 
the full second order wave equation or Kirchoff momentum equation. 

a transmission line based investigation and validation of MacCormack's 
characteristic factorization algorithm. 

an experimental program to measure information propagation effects in 
transmission lines and Stirling machines using precise digital acoustic 
techniques. These tests should emphasize the accurate measurement of 
phase angles in order to validate this aspect of Iberall's analysis as 
well as to determine the information propagation characteristics of 
actual Stirling hardware. 

The results of part B suggest that the false diffusion correction 
methodology developed may yield some improvement in the accuracy of existing 
one -dimensional Stirling machine simulation codes. This improvement would 
result from standard steady-state heat transfer correlations yielding better 
accuracy while the enthalpy fluxes at the heater and cooler entrance and exit 
ports also would be more accurately modelled. 

Inclusion of entropy transport (as opposed to entropy generation alone) 
in Stirling machine codes would yield a reliable method for tracking 
irreversibilities within the engine and would allow a better comparison of the 
performance of different simulation codes. Hence, as a next step, it seems 
appropriate that the false diffusion correction algorithm and entropy 
transport equation be ported to existing codes which should then be re- 
validated against experimental data. 
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